#########################
#### Call libraries. ####

library(ssanv)
library(exactci)
library(exact2x2)
library(sensitivitymw)

library(iterators)
library(foreach)
library(doParallel)
library(doMC)

library(plyr)


###################################
#### Read in the matched data. ####

trad_gensurg <- read.csv("/Users/sharpe/Desktop/Young Surgeons/data/simul_trad_patient_match_gensurg_11.2.2018.csv")
mod_gensurg <- read.csv("/Users/sharpe/Desktop/Young Surgeons/data/simul_mod_patient_match_gensurg_11.2.2018.csv")
trad_ortho <- read.csv("/Users/sharpe/Desktop/Young Surgeons/data/simul_trad_patient_match_ortho_11.2.2018_exact_procs.csv")
mod_ortho <- read.csv("/Users/sharpe/Desktop/Young Surgeons/data/simul_mod_patient_match_ortho_11.2.2018_exact_procs.csv")

trad_gensurg <- trad_gensurg[,-which(names(trad_gensurg) %in% c("total_cost_complete", "total_cost_index_only","icu_days","total_cost_other"))]
mod_gensurg <- mod_gensurg[,-which(names(mod_gensurg)%in% c("total_cost_complete", "total_cost_index_only","icu_days","total_cost_other"))]
trad_ortho <- trad_ortho[,-which(names(trad_ortho)%in% c("total_cost_complete", "total_cost_index_only","icu_days","total_cost_other"))]
mod_ortho <- mod_ortho[,-which(names(mod_ortho)%in% c("total_cost_complete", "total_cost_index_only","icu_days","total_cost_other"))]


###############################################
#### Read in the updated icu_yes variable. ####

trad_ortho_icu <- read.csv("/Users/sharpe/Desktop/Young Surgeons/data/trad_ortho_new_cost_and_icu.csv")
mod_ortho_icu <- read.csv("/Users/sharpe/Desktop/Young Surgeons/data/mod_ortho_new_cost_and_icu.csv")
trad_gensurg_icu <- read.csv("/Users/sharpe/Desktop/Young Surgeons/data/trad_gensurg_new_cost_and_icu.csv")
mod_gensurg_icu <- read.csv("/Users/sharpe/Desktop/Young Surgeons/data/mod_gensurg_new_cost_and_icu.csv")

trad_ortho_icu <- trad_ortho_icu[,c("newID_1","newID_10","total_cost_complete", "total_cost_index_only","icu_days","total_cost_other")]
mod_ortho_icu <- mod_ortho_icu[,c("newID_1","newID_10","total_cost_complete", "total_cost_index_only","icu_days","total_cost_other")]
trad_gensurg_icu <- trad_gensurg_icu[,c("newID_1","newID_10","total_cost_complete", "total_cost_index_only","icu_days","total_cost_other")]
mod_gensurg_icu <- mod_gensurg_icu[,c("newID_1","newID_10","total_cost_complete", "total_cost_index_only","icu_days","total_cost_other")]

#####################################
#### Merge the new ICU variable. ####

trad_gensurg <- merge(x=trad_gensurg, y=trad_gensurg_icu, by=c("newID_1","newID_10"), all.x=T)
mod_gensurg <- merge(x=mod_gensurg, y=mod_gensurg_icu, by=c("newID_1","newID_10"), all.x=T)
trad_ortho <- merge(x=trad_ortho, y=trad_ortho_icu, by=c("newID_1","newID_10"), all.x=T)
mod_ortho <- merge(x=mod_ortho, y=mod_ortho_icu, by=c("newID_1","newID_10"), all.x=T)

gensurg_match <- rbind.fill(trad_gensurg, mod_gensurg)
ortho_match <- rbind.fill(trad_ortho, mod_ortho)

gensurg_match$pair_id_pats <- paste(gensurg_match$pair_id_pats,"_g",gensurg_match$era_trad,sep="")
ortho_match$pair_id_pats <- paste(ortho_match$pair_id_pats,"_o",ortho_match$era_trad, sep="")


both_match <- rbind.fill(gensurg_match, ortho_match)

trad_match <- rbind.fill(gensurg_match[which(gensurg_match$era_trad==1),], ortho_match[which(ortho_match$era_trad==1),])
mod_match <- rbind.fill(gensurg_match[which(gensurg_match$era_mod==1),], ortho_match[which(ortho_match$era_mod==1),])

trad_gensurg <- subset(gensurg_match, era_trad==1)
trad_ortho <- subset(ortho_match, era_trad==1)
mod_gensurg <- subset(gensurg_match, era_mod==1)
mod_ortho <- subset(ortho_match, era_mod==1)

#####################################
#### Read in the unmatched data. ####

mod_ortho_unmatched <- read.csv("/Users/sharpe/Desktop/Young Surgeons/data/transfer_mod_ortho_simul_11_02.csv")
trad_ortho_unmatched <- read.csv("/Users/sharpe/Desktop/Young Surgeons/data/transfer_trad_ortho_simul_11_02.csv")
mod_gensurg_unmatched <- read.csv("/Users/sharpe/Desktop/Young Surgeons/data/transfer_mod_gensurg_simul_11_02.csv")
trad_gensurg_unmatched <- read.csv("/Users/sharpe/Desktop/Young Surgeons/data/transfer_trad_gensurg_simul_11_02.csv")

trad_gensurg_unmatched <- trad_gensurg_unmatched[,-which(names(trad_gensurg_unmatched) %in% c("total_cost_complete", "total_cost_index_only","icu_days","total_cost_other"))]
mod_gensurg_unmatched <- mod_gensurg_unmatched[,-which(names(mod_gensurg_unmatched)%in% c("total_cost_complete", "total_cost_index_only","icu_days","total_cost_other"))]
trad_ortho_unmatched <- trad_ortho_unmatched[,-which(names(trad_ortho_unmatched)%in% c("total_cost_complete", "total_cost_index_only","icu_days","total_cost_other"))]
mod_ortho_unmatched <- mod_ortho_unmatched[,-which(names(mod_ortho_unmatched)%in% c("total_cost_complete", "total_cost_index_only","icu_days","total_cost_other"))]

trad_gensurg_unmatched <- merge(x=trad_gensurg_unmatched, y=trad_gensurg_icu, by=c("newID_1","newID_10"), all.x=T)
mod_gensurg_unmatched <- merge(x=mod_gensurg_unmatched, y=mod_gensurg_icu, by=c("newID_1","newID_10"), all.x=T)
trad_ortho_unmatched <- merge(x=trad_ortho_unmatched, y=trad_ortho_icu, by=c("newID_1","newID_10"), all.x=T)
mod_ortho_unmatched <- merge(x=mod_ortho_unmatched, y=mod_ortho_icu, by=c("newID_1","newID_10"), all.x=T)


mod_ortho_unmatched$new_surgeon_flag <- ifelse(mod_ortho_unmatched$Mod_New_Surgeon_flag_numeric==1,1,0)
trad_ortho_unmatched$new_surgeon_flag <- ifelse(trad_ortho_unmatched$Trad_New_Surgeon_flag_numeric==1,1,0)
mod_gensurg_unmatched$new_surgeon_flag <- ifelse(mod_gensurg_unmatched$Mod_New_Surgeon_flag_numeric==1,1,0)
trad_gensurg_unmatched$new_surgeon_flag <- ifelse(trad_gensurg_unmatched$Trad_New_Surgeon_flag_numeric==1,1,0)

ortho_unmatched <- rbind(mod_ortho_unmatched, trad_ortho_unmatched)
gensurg_unmatched <- rbind(mod_gensurg_unmatched, trad_gensurg_unmatched)

trad_unmatched <- rbind.fill(trad_ortho_unmatched, trad_gensurg_unmatched)
mod_unmatched <- rbind.fill(mod_ortho_unmatched, mod_gensurg_unmatched)

all_unmatched <- rbind.fill(mod_ortho_unmatched, mod_gensurg_unmatched, trad_ortho_unmatched, trad_gensurg_unmatched)

#########################################################
#### Perform M-estimates on the continuous outcomes. ####


#### Define the continuous outcome variables we want to calculate M-estimates for. ####

con_outcomes <- c("total_cost_complete","total_cost_index_only","anesth_time_claim","icu_days","los_long", "basic_los")

#### Create an M-estimate function for the matched pairs. ####

mestimates <- function(data, treatment, pair_id, variable, comparison){
  
  test_data_t1 <- data[which(data[,treatment]==1),c(pair_id,variable)]
  test_data_t0 <- data[which(data[,treatment]==0),c(pair_id,variable)]
  
  colnames(test_data_t1)[2] <- c("var_t1")
  colnames(test_data_t0)[2] <- c("var_t0")
  
  test_data <- merge(x=test_data_t1, y=test_data_t0, by=c(pair_id))
  test_data <- as.matrix(subset(test_data, select=c("var_t1","var_t0")))
  
  m1 <- senmwCI(y=test_data_t1[,c("var_t1")], gamma=1, method="h", inner=0, trim=1, lambda=0.98, m=1, m1=1, m2=1, alpha=0.05, one.sided=F, tol=NULL, interval=NULL, detail=T)
  m0 <- senmwCI(y=test_data_t0[,c("var_t0")], gamma=1, method="h", inner=0, trim=1, lambda=0.98, m=1, m1=1, m2=1, alpha=0.05, one.sided=F, tol=NULL, interval=NULL, detail=T)
  diff <- senmwCI(y=(test_data_t1[,c("var_t1")] - test_data_t0[,c("var_t0")]), gamma=1, method="h", inner=0, trim=1, lambda=0.98, m=1, m1=1, m2=1, alpha=0.05, one.sided=F, tol=NULL, interval=NULL, detail=T)
  
  if(diff$PointEstimate[1] < 0){
    
    test_data <- test_data*(-1)
    
  }
  
  
  m <- senmw(y=test_data, gamma=1, method="h", inner=0, trim=1, lambda=0.98, tau=0, m=1, m1=1, m2=1)
  

  
  m_vec <- c(comparison, variable, (m$pval)*2, m$deviate, m$statistic,m$expectation, m$variance, m1$PointEstimate[1], m1$Confidence.Interval[1],m1$Confidence.Interval[2],m0$PointEstimate[1], m0$Confidence.Interval[1], m0$Confidence.Interval[2], diff$PointEstimate[1], diff$Confidence.Interval[1], diff$Confidence.Interval[2])
  names(m_vec) <- c("Comparison","Outcome","M Statistic P-value", "Deviate", "Statistic", "Expectation", "Variance", "New Surgeons M-estimate", "New Surgeons Lower 95% CI", "New Surgeons Upper 95% CI", "Exp. Surgeons M-estimate", "Exp. Surgeons Lower 95% CI", "Exp. Surgeons Upper 95% CI", "New - Exp. M-estimate",
                    "New - Exp. Lower 95% CI", "New - Exp. Upper 95% CI")
  
  
  
  return(m_vec)
  
  
}

#### Create an M-estimates function for overall estimates. ####

mestimates_overall <- function(data, treatment, variable, comparison){
  
  test_data_t0 <- data[which(data[,treatment]==0),c(variable)]
  test_data_t1 <- data[which(data[,treatment]==1),c(variable)]
  

  
  m0 <- senmwCI(y=test_data_t0, gamma=1, method="h", inner=0, trim=1, lambda=0.98, m=1, m1=1, m2=1, alpha=0.05, one.sided=F, tol=NULL, interval=NULL, detail=T)
  m1 <- senmwCI(y=test_data_t1, gamma=1,method="h",inner=0,trim=1,lambda=0.98,m=1,m1=1,m2=1,alpha=0.05,one.sided=F,tol=NULL, interval=NULL, detail=T)
  m_vec <- c(comparison, variable, m0$PointEstimate[1], m0$Confidence.Interval[1],m0$Confidence.Interval[2],  m1$PointEstimate[1], m1$Confidence.Interval[1],m1$Confidence.Interval[2])
  names(m_vec) <- c("Comparison","Outcome","Exp. Surgeon M-Estimate", "Exp Surgeon - Lower 95% CI", "Exp Surgeon - Upper 95% CI","New Surgeon M-Estimate", "New Surgeon - Lower 95% CI", "New Surgeon - Upper 95% CI")
  
  return(m_vec)
}


#### Loop through the datasets and the variables and save the output. ####


# mest_both <- NULL
# mest_trad <- NULL
# mest_mod <- NULL
# mest_gensurg <- NULL
# mest_ortho <- NULL
# mest_trad_ortho <- NULL
# mest_trad_gensurg <- NULL
# mest_mod_ortho <- NULL
# mest_mod_gensurg <- NULL

mest_vec_all <- NULL

matched_list <- vector(length=9, mode="list")
matched_list[[1]] <- both_match
matched_list[[2]] <- trad_match
matched_list[[3]] <- mod_match
matched_list[[4]] <- gensurg_match
matched_list[[5]] <- ortho_match
matched_list[[6]] <- trad_ortho
matched_list[[7]] <- trad_gensurg
matched_list[[8]] <- mod_ortho
matched_list[[9]] <- mod_gensurg

start <- Sys.time()
registerDoMC(cores=10)



mest_vec_all <- foreach(i=1:length(con_outcomes),.combine='rbind') %:%
                  foreach(j=1:length(matched_list),.combine='rbind') %dopar% {
                    
                    comparison <- ifelse(j==1,"all matched data",ifelse(j==2,"trad",ifelse(j==3,"mod",ifelse(j==4,"gensurg",ifelse(j==5,"ortho",ifelse(j==6,"trad ortho",ifelse(j==7,"trad gensurg",ifelse(j==8,"mod ortho","mod gensurg"))))))))
                    
                    a_vec <- mestimates(data=matched_list[[j]], treatment="new_surgeon_flag", pair_id="pair_id_pats", variable=con_outcomes[i], comparison=comparison)
                    
                    
                    
}

end <- Sys.time()
total_time1 <- end-start
total_time1

write.csv(x=mest_vec_all, file="/Users/sharpe/Desktop/Young Surgeons/table/simul_matched_m_estimates_REDONE_11.28.2018.csv")


mest_vec_overall <- NULL

unmatched_list <- vector(length=5,mode="list")
unmatched_list[[1]] <- gensurg_unmatched
unmatched_list[[2]] <- ortho_unmatched
unmatched_list[[3]] <- trad_unmatched
unmatched_list[[4]] <- mod_unmatched
unmatched_list[[5]] <- all_unmatched



mest_vec_overall <- foreach(i=1:length(con_outcomes),.combine='rbind') %:%
                      foreach(j=1:length(unmatched_list),.combine='rbind') %dopar% {
                        
                        comparison=ifelse(j==1,"all gensurg",ifelse(j==2,"all ortho",ifelse(j==3,"all trad",ifelse(j==4,"all mod","all"))))
                        b_vec <- mestimates_overall(data=unmatched_list[[j]], treatment="new_surgeon_flag",variable=con_outcomes[i], comparison=comparison)
                        
}

write.csv(x=mest_vec_overall, file="/Users/sharpe/Desktop/Young Surgeons/table/simul_m_estimates_overall_REDONE_11.28.2018.csv")




